########### TABLE OF CONTENTS ############################
# 0. Libraries - define them all at the top
# 1. Define Functions
# 2. Human inputs - manually enter file name as variable f (USER PROVIDES file name)
# 3. Populate alldata variable from file
# 4a. Centre, Filter/Smooth the x,y detected data
# 4b. Simulate x,y data to check calculations - Remove for actual analysis
# 5. Clean up the Fish Side Detection
# 6. Define temperature fish selects and subset by side
# 7. Movement trajectories using bearing function
# 8. Plot Fish smooth (x,y) dataa
# 9. Plot the trajectory data, polar plot, histograms
#10. Summaries and final calculations
########### End of TOC ##################################
#0. Libraries required ####
# if missing, type install.packages("packagename") in console
library(plotrix)
library(Thermimage)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(signal)
##
## Attaching package: 'signal'
## The following objects are masked from 'package:stats':
##
## filter, poly
library(pastecs)
## Loading required package: boot
library(graphics)
library(fields)
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.1-1 (2017-07-02) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: maps
##
## Attaching package: 'fields'
## The following object is masked from 'package:plotrix':
##
## color.scale
library(plyr)
##
## Attaching package: 'plyr'
## The following object is masked from 'package:maps':
##
## ozone
source("https://www.dropbox.com/s/168hsu4xuus0o1d/CommonFunctions.R?dl=1")
## Loading required package: Matrix
## Loading required package: MASS
##
## arm (Version 1.9-3, built: 2016-11-21)
## Working directory is /Users/GlennTattersall/Dropbox/R/MyProjects/ActivityDirection
##
## Attaching package: 'arm'
## The following object is masked from 'package:spam':
##
## display
## The following object is masked from 'package:boot':
##
## logit
## The following object is masked from 'package:plotrix':
##
## rescale
bias.index.c<-NULL
options(warn=-1)
#1. Define Functions ####
bearing<-function(x,y,zmin=0)
# this function accepts x and y values of detected coordinates for movement
# and returns a data frame of four variables: magnitude, theta, change, alpha
# magnitude and theta define the vectors describing the direction of movement
# between time points i and i+1. By necessity, magnitude and theta will be
# 1 row shorter than the x,y input coordinates
# change and alpha represent the change in magnitude and direction in the vector
# of movement
# zmin is the lowest acceptable level of magnitude required to accept the vector
# of movement of being realisitc (i.e. zmin is the minimum detectable magnitude)
# http://stackoverflow.com/questions/2571160/calculate-direction-angle-from-two-vectors
{
dx<-diff(x)
dy<-diff(y)
velocity<-sqrt(dy^2+dx^2)
# remove really small changes (i.e. zmin=3 units or less)
dx[which(abs(velocity)<zmin)]<-NA
dy[which(abs(velocity)<zmin)]<-NA
velocity[which(abs(velocity)<zmin)]<-NA
theta<-180*atan2(dy,dx)/pi
# atan2 will calculate angles with respect to the horizontal.
alpha<-diff(theta)
alpha[which(alpha>180)]<-(-360+alpha[which(alpha>180)])
alpha[which(alpha<(-180))]<-360+alpha[which(alpha<(-180))]
# alpha is the instantaneous change in theta with respect to index
# alternatively called the relative angle or turning angle
# if the index is in unitary measures of time, alpha becomes degrees/time
# a +ve alpha corresponds to a left turn, or counter-clockwise turn
# a -ve alpha corresponds to a right turn, or clock-wise turn
acceleration<-diff(velocity)
# change is the magnitude of the relative turn
alpha<-c(NA,alpha)
acceleration<-c(NA,acceleration)
# assign NA to the first value as a placeholder, so that alpha and change
# align with the n+1 data points
result<-data.frame(velocity, theta, acceleration, alpha)
}
addgradient <- function()
{
# creates the gradient fill colour ramp for smoothScatter plots
xm <- get('xm', envir = parent.frame(1))
ym <- get('ym', envir = parent.frame(1))
z <- get('dens', envir = parent.frame(1))
colramp <- get('colramp', parent.frame(1))
image.plot(xm,ym,z, col = colramp(256), legend.only = T, add =F)
}
spiral<-function(t=seq(1,100,0.1), size=5, direction="left")
{
if(direction=="right")
{
x<-size*t*sin(t)
y<-size*t*cos(t)
}
else if(direction=="left")
{
x<-size*t*cos(t)
y<-size*t*sin(t)
}
out<-cbind(x,y)
}
alphacumulate<-function(alphas, events){
alphas[is.na(alphas)]<-0
cumulate<-NULL
z<-max(events)
for(i in 1:z){
tmp<-cumsum(alphas[events==i])
cumulate<-c(cumulate, tmp)
}
return(cumulate)
}
#2. Human Inputs - Filename ####
mainDir<-"~/Dropbox/R/MyProjects/ActivityDirection/data"
outputDir<-"~/Dropbox/R/MyProjects/ActivityDirection/output"
#mainDir<-"~/Desktop/Justin/data/" # Folder for student
#outputDir<-"~/Desktop/Justin/output/" # Folder for student
dir.create(outputDir, showWarnings = F, recursive = FALSE, mode = "0777")
setwd(mainDir)
l.files<-list.files()
l.files
## [1] "22518 - 19SEB19.csv" "22818 - 18SEB11.csv" "3118 - 15LSJ11.csv"
# User input ####
f<-l.files[3]
f
## [1] "3118 - 15LSJ11.csv"
#for(f in l.files)
#{
setwd(mainDir)
Rdatafileout<-paste(sub(".csv", "", f),".rds", sep="")
timesummaryout<-paste(sub(".csv", "", f),"_by_time.csv", sep="")
shuttlesummaryout<-paste(sub(".csv", "", f),"_by_shuttle.csv", sep="")
fileout<-"Output.csv"
#3. Populate alldata data.frame ####
alldata<-read.csv(f, header=TRUE, skip=0, row.names=NULL)
alldata<-alldata[c("Real", "Elapsed", "X.Pos", "Y.Pos", "Area..px.", "Temp..Left.",
"Temp..Right.", "Fish.Side", "Distance.m.", "Distance.p.",
"Velocity.m.s.", "Velocity.p.s.")]
#alldata<-alldata[,-c(16:22)] # remove columns added by Nicole ####
alldata<-na.omit(alldata)
colnames(alldata)<-c("Real", "Elapsed", "X.Pos", "Y.Pos", "Area", "Temp.Left",
"Temp.Right", "Fish.Side", "Distance.m", "Distance.p",
"Velocity.m.s", "Velocity.p.s")
# Look at Area size. If Area=0, problem with detection
ind<-which(alldata$Area==0)
# replace the erroneous detections with the previous row.
alldata$X.Pos[ind]<-alldata$X.Pos[(ind-1)]
alldata$Y.Pos[ind]<-alldata$Y.Pos[(ind-1)]
alldata$X.Pos[which(alldata$X.Pos<0)]<-NA
alldata$Y.Pos[which(alldata$Y.Pos<0)]<-NA
alldata$X.Pos<-na.fill(alldata$X.Pos, "extend")
alldata$Y.Pos<-na.fill(alldata$Y.Pos, "extend")
# par(mfrow=c(1,1), mar=c(10,5,5,5))
# plot(alldata$Temp.Left, main="Click at the time where ramping starts, then <ESC>")
# title(sub="Only the data following this time point will be analysed", line=4)
oldpar<-par()
#
# timeremove<-0
# timeremove<-round(mean(locator(1)$x),0)
# # rem out line above with timeremove and locator if running Knit #####
#
# if(timeremove<1) timeremove<-1
# cat("The first", timeremove, "data points will be removed")
# cat("\n")
# # User input #####
# alldata<-alldata[-c(1:timeremove),]
# # remove the first hour (i.e. 3600 s) of data.
#4a. Centre, Filter, and Smooth the x and y data ####
# alldata$X.Pos<-alldata$X.Pos-(max(alldata$X.Pos,na.rm=TRUE) +
# (min(alldata$X.Pos,na.rm=TRUE))/2
# alldata$Y.Pos<-alldata$Y.Pos-(max(alldata$Y.Pos,na.rm=TRUE) +
# min(alldata$Y.Pos,na.rm=TRUE))/2
alldata$X.Pos<-alldata$X.Pos-160
alldata$Y.Pos<-alldata$Y.Pos-120
half.x<-0
# centre the x data, assuming that the max and min values detected are the extremes
# Convert the X.Pos into true distance (cm)
fit<-lm(100*Distance.m ~ Distance.p, data=alldata)
b<-fit$coef[1]
b<-0
m<-fit$coef[2]
# use the data file to ascertain correlation between pixel and cm
alldata$X.Pos<-b+m*alldata$X.Pos
alldata$Y.Pos<-b+m*alldata$Y.Pos
alldata<-alldata[c("Real", "Elapsed", "X.Pos", "Y.Pos", "Area", "Temp.Left",
"Temp.Right", "Fish.Side")]
# convert X.Pos and Y.Pos into actual distance units (cm)
alldata$Real<-as.POSIXct(as.character(alldata$Real), tz="", "%H:%M:%S")
alldata$Elapsed<-as.POSIXct(as.character(alldata$Elapsed), tz="", "%H:%M:%S")
Time.diff<-as.integer(alldata$Elapsed-alldata$Elapsed[1])
t<-1:nrow(alldata)
alldata<-data.frame(Time.diff, alldata, row.names=NULL)
x<-alldata$X.Pos
y<-alldata$Y.Pos
bf<-butter(2, 0.2, type="low")
x.filt1<-filter(bf, x)
y.filt1<-filter(bf, y)
x.filt2<-sgolayfilt(x, 9, 55)
y.filt2<-sgolayfilt(y, 9, 55)
x.filt<-as.numeric((x.filt1+x.filt2)/2)
y.filt<-as.numeric((y.filt1+y.filt2)/2)
# User input ####
# empirically determined filters - may require adjustment
x.filt<-na.locf(x.filt, na.rm=FALSE)
y.filt<-na.locf(y.filt, na.rm=FALSE)
# replace any NA values with previous non-NA value
par(mfrow=c(1,2))
plot(x[10:100], type="l", main="Sample of x values",
xlab="Time", ylab="x", col="grey", lty=2)
lines(x.filt[10:100], type="l", col="blue")
legend(x="topright", c("Raw", "Filtered"), col=c("grey", "blue"), lty=c(2,1), bty="n")
plot(y[10:100], type="l", main="Sample of y values",
xlab="Time", ylab="y", col="grey", lty=2)
lines(y.filt[10:100], type="l", col="blue")
legend(x="topright", c("Raw", "Filtered"), col=c("grey", "blue"), lty=c(2,1), bty="n")
# 10 data points shown for illustrative purposes")
alldata<-data.frame(alldata, x.filt, y.filt, row.names=NULL)
plot(x, y, main="Raw Data",
type='l', xlim=c(min(x),max(x)), ylim=c(min(y),max(y)))
plot(x.filt, y.filt, main="Filtered Data",
type='l', xlim=c(min(x),max(x)), ylim=c(min(y),max(y)))
#4b. Simulate x,y as a spiral or random #####
rand<-2
# set rand to 1 to run a random sample, #####
# set rand to 0 to run a spiral,
# set rand to any other number to use imported data
maxrepeats<-1
if(rand==1)
{
maxrepeats<-10000
}
#for(m in 1:maxrepeats)
#{
if(rand==1)
{
# theta.rand<-runif(nrow(alldata), -1, 1)
# mag.rand<-runif(nrow(alldata), 0, 1)
# x.rand<-mag.rand*sin(theta.rand)
# x.rand<-cumsum(x.rand)
# y.rand<-mag.rand*cos(theta.rand)
# y.rand<-cumsum(y.rand)
#
index.rand1<-sample(1:nrow(alldata), nrow(alldata))
#index.rand2<-sample(1:nrow(alldata), nrow(alldata))
x.rand<-alldata$x.filt[index.rand1]
y.rand<-alldata$y.filt[index.rand1]
alldata$x.filt<-x.rand
alldata$y.filt<-y.rand
}
if(rand==0)
{
spir<-spiral(seq(1,nrow(alldata)), 1, "right")
alldata$x.filt<-spir[,1]
alldata$y.filt<-spir[,2]
}
#5. Clean up Fish.Side ??? detections ####
alldata$Fish.Side<-as.character(alldata$Fish.Side)
q<-which(!alldata$Fish.Side=="???")
if(alldata$Fish.Side[1]=="???") alldata$Fish.Side[1]<-alldata$Fish.Side[q[1]]
for(i in 1:nrow(alldata))
{
j<-i-1
#if(alldata$Fish.Side[i]=="???")
# {
if(x.filt[i]<half.x) alldata$Fish.Side[i]<-"LEFT"
if(x.filt[i]>half.x) alldata$Fish.Side[i]<-"RIGHT"
if(x.filt[i]==half.x) alldata$Fish.Side[i]<-alldata$Fish.Side[j]
# }
}
alldata$Fish.Side<-as.factor(alldata$Fish.Side)
#6. Define fish temperature ####
slopeleft<-slopeEveryN(alldata$Temp.Left, 30)
sloperight<-slopeEveryN(alldata$Temp.Right, 30)
overallslope<-(slopeleft+sloperight)/2
left.side<-which(alldata$Fish.Side=="LEFT")
right.side<-which(alldata$Fish.Side=="RIGHT")
fish.temp<-1:nrow(alldata) # populate variable fish.temp
fish.temp[left.side]<-alldata$Temp.Left[left.side]
fish.temp[right.side]<-alldata$Temp.Right[right.side]
alldata<-data.frame(alldata, fish.temp, row.names=NULL)
fs<-as.integer(alldata$Fish.Side)-1
# convert Fish.Side character to an integer vector, where Left=0, Right=1
sideswitch<-cbind(1:nrow(alldata), c(0,diff(fs)))
sideswitch<-data.frame(sideswitch)
colnames(sideswitch)<-c("Index","Switch")
# create a indexed matrix that indicates where shuttles occurred
# +1 = moving from left into right (i.e. lower escape)
# - = moving from right into left (i.e. upper escape)
sideswitch<-subset(sideswitch, Switch >0 | Switch<0)
sideswitch$Event<-1:nrow(sideswitch)
le.index<-sideswitch[which(sideswitch[,2]==1), 1]
ue.index<-sideswitch[which(sideswitch[,2]==-1), 1]
let<-alldata$Temp.Left[le.index]
uet<-alldata$Temp.Right[ue.index]
leshutNo<-1:length(le.index)
ueshutNo<-1:length(ue.index)
alldata$Event<-NA
alldata$Event[sideswitch$Index]<-sideswitch$Event
if(is.na(alldata$Event[1])){
alldata$Event<-alldata$Event+1
alldata$Event[1]<-1
}
alldata$Event<-na.locf(alldata$Event)
#7. Work out movement trajectories using bearing function #####
fish.bearing<-bearing(alldata$x.filt, alldata$y.filt, 0)
#fish.bearing<-bearing(alldata$X.Pos,alldata$Y.Pos, 0)
delta.time<-diff(alldata$Time.diff)
velocity<-fish.bearing[,1]/delta.time
# divide by delta.time in case each row is not 1 sec apart
velocity[is.infinite(velocity)]<-NA
# replace INF values with NA
velocity<-na.locf(velocity, na.rm=FALSE)
# use na.locf to replace NA values with previous entry
theta<-fish.bearing[,2]/delta.time
# divide by alldata$Time.diff in case each row is not 1 sec apart
theta[is.infinite(theta)]<-NA
theta<-na.locf(theta, na.rm=FALSE)
# use na.locf to replace NA values with previous entry
# this will bestow memory to the direction the fish was last moving
acceleration<-abs(fish.bearing[,3]/delta.time)
acceleration[is.infinite(acceleration)]<-NA
alpha<-fish.bearing[,4]/delta.time
alpha[is.infinite(alpha)]<-NA
# divide by alldata$Time.diff in case each row is not 1 sec apart
acceleration<-na.locf(acceleration, na.rm=FALSE)
alpha<-na.locf(alpha, na.rm=FALSE)
# alphas = 0 will mean a stationary animal
# alphas <> 0 can then be examined for how large or directional they are
alpha[which(alpha>180)]<-(-360+alpha[which(alpha>180)])
alpha[which(alpha<(-180))]<-360+alpha[which(alpha<(-180))]
# convert alpha angles so that 0 to 180 = left, counterclockwise
# and -0 to -180 = right, clockwise
#acceleration<-c(NA, acceleration)
#alpha<-c(NA, alpha)
alldata<-data.frame(alldata[-1,], delta.time, velocity, theta, acceleration, alpha, row.names=NULL)
# remove first row from original dataframe before adding the bearing output
# since the bearing data is based on one less data point
colnames(alldata)<-c("Time.diff", "Real", "Elapsed", "X.Pos", "Y.Pos", "Area", "Temp.Left",
"Temp.Right", "Fish.Side", "x.filt", "y.filt", "fish.temp", "Event",
"Delta.time", "velocity", "theta", "acceleration", "alpha")
alldata$distance<-alldata$velocity*alldata$Delta.time
localmax<-function(dat, ind, rng=9){
startind<-ind-(rng-1)/2
endind<-ind+(rng-1)/2
allindex<-unlist(Map(':', startind, endind))
z<-matrix(dat[allindex], nrow=rng)
lmax<-apply(z, 2, max)
return(lmax) # returns the local max value
}
# Extract velocity and acceleration values coinciding with each shuttle. extract the maximum value to
# estimate burst values and to avoid stationary data points
vellet<-localmax(alldata$velocity, ind=le.index, rng=3)
veluet<-localmax(alldata$velocity, ind=ue.index, rng=3)
acclet<-localmax(alldata$acceleration, ind=le.index, rng=3)
accuet<-localmax(alldata$acceleration, ind=ue.index, rng=3)
# create shuttle only data frames for eventual export
leshuttles<-data.frame(1:length(let), le.index, alldata$Time.diff[le.index], let, vellet, acclet)
colnames(leshuttles)<-c("Shuttle.No", "LE.Index", "LE.Time", "LET", "Veloc.LET", "Accel.LET")
ueshuttles<-data.frame(1:length(uet), ue.index, alldata$Time.diff[ue.index], uet, veluet, accuet)
colnames(ueshuttles)<-c("Shuttle.No", "UE.Index", "UE.Time", "UET", "Veloc.UET", "Accel.UET")
allshuttles<-merge(leshuttles, ueshuttles)
#plot(allshuttles$Veloc.LET~allshuttles$LET)
#plot(allshuttles$Veloc.UET~allshuttles$UET)
#8. Plot Fish smoothed (X,Y) data ####
par(mfrow=c(1,1), mar=c(5,5,5,5))
plot(alldata$X.Pos, alldata$Y.Pos, type="l",
main="Fish Positions - Raw", xlab="x", ylab="y")
plot(alldata$x.filt,alldata$y.filt, type="l",
main="Fish Positions - Smoothed", xlab="x", ylab="y")
#smooth.palette <- colorRampPalette(blues9[-(1:2)])
smooth.palette <- palette.choose("glowbow")
smoothScatter(alldata$X.Pos, alldata$Y.Pos, nbin=256, bty="n", xlab="X Position (cm)",
ylab="Y Position (cm)",
main="Density Plot of Fish Positions - Raw", colramp=heat.colors,
xlim=c(min(alldata$X.Pos),max(alldata$X.Pos)),
ylim=c(min(alldata$Y.Pos),max(alldata$Y.Pos)),
postPlotHook=addgradient)
smoothScatter(alldata$x.filt, alldata$y.filt, nbin=256, bty="n",xlab="X Position (cm)",
ylab="Y Position (cm)",
main="Density Plot of Fish Positions - Smoothed", colramp=heat.colors,
xlim=c(min(alldata$X.Pos),max(alldata$X.Pos)),
ylim=c(min(alldata$Y.Pos),max(alldata$Y.Pos)),
postPlotHook=addgradient)
# Plot the X,Y positions, colour coded by alpha
par(oldpar)
par(mfrow=c(1,1), mar=c(5,5,5,5))
plot(0, 0, type="n", pch=16, xlim=c(min(alldata$x.filt),1.1*max(alldata$x.filt)),
ylim=c(min(y.filt),1.1*max(y.filt)), main="Detected Positions",
xlab="X", ylab="Y")
alphaPal <- colorRampPalette(c('red', 'white', 'blue'))
a.col<-alldata$alpha
a.col[which(a.col>30)]<-30
a.col[which(a.col<(-30))]<--30
alphaCol <- alphaPal(10)[as.numeric(cut(a.col, breaks = 10))]
# points(alldata$x.filt[-1], alldata$y.filt[-1], cex=0.5, pch = 21, col=alphaCol,
# bg=alphaCol)
points(alldata$X.Pos[-1], alldata$Y.Pos[-1], cex=0.5, pch = 21, col=alphaCol,
bg=alphaCol)
legend(x="topright", c("Right turn", "Left turn"), col=c("red", "blue"), pch=c(1,1),
bty="n")
#9. Plot the trajectory data ####
par(mfrow=c(1,2))
polar.plot(velocity, theta, line.col="red", point.col="red", rp.type="s",
start=0, main=("Movement (° vs. cm/s) of Fish"), cex=0.2)
title(sub="Angle (°) is relative to horizontal = 0°", line=1)
polar.plot(acceleration[-1], alpha[-1], point.col="red", start=90, rp.type="s",
main="Turning Angle (° vs. cm/s/s)", cex=0.2,
labels=c(0,20,40,60,80,100,120,140,160,180,-160,-140,-120,-100,
-80,-60,-40,-20))
title(sub="+° = left, -° = right", line=1)
bin.labels<-seq(-172.5,172.5,15)
bins<-seq(-180,180,15)
par(mfrow=c(1,1))
freq.alpha<-as.numeric(table(cut(alpha, breaks=bins, labels=bin.labels)))
polar.plot(100*freq.alpha/sum(freq.alpha), bin.labels, line.col="black",
rp.type="r", start=90, lwd=10, lend=1,
main="Distribution of Alphas (% of observed turning angles)",
labels=c(0,20,40,60,80,100,120,140,160,180,-160,-140,-120,-100,
-80,-60,-40,-20),
mar=c(5,5,5,5))
title(sub="+° = left, -° = right", line=1)
par(oldpar)
par(mfrow=c(2,2), mar=c(5,5,2,5))
hist(velocity, main="Velocity", xlab="Velocity (cm/s)")
hist(theta, main="Thetas", xlab="Angle (°)", ylab="# of observations")
hist(acceleration[-1], main="Acceleration", xlab="Acceleration (cm/s/s)",
ylab="# of observations")
hist(alpha[-1], main="Alphas", xlab="Angle (°)", ylab="# of observations" )
mean(alpha, na.rm=TRUE)
## [1] -0.9023
par(mfrow=c(1,2), mar=c(5,5,5,5))
hist(na.omit(alpha), main="Frequency of all alphas", xlab="Turning angle (°)",
ylab="# of observations")
hist(na.omit(alpha[which(abs(alpha)>20)]), main="Frequency of | alphas | > 20°",
xlab="Turning angle (°)", ylab="# of observations")
par(mfrow=c(1,1))
freq.alpha20<-as.numeric(table(cut(alpha[which(abs(alpha)>20)], breaks=bins,
labels=bin.labels)))
polar.plot(100*freq.alpha20/sum(freq.alpha20), bin.labels, line.col="black",
rp.type="r", start=90, lwd=10, lend=1,
main="Distribution of | alphas | > 20° (% of observed turning angles)",
labels=c(0,20,40,60,80,100,120,140,160,180,-160,-140,-120,-100,
-80,-60,-40,-20),
mar=c(5,5,5,5))
title(sub="+° = left, -° = right", line=1)
# Cumulative Alpha calculation ####
# calculative cumulative alpha and reset every time the fish switches sides (event)
alldata$cumalpha<-alphacumulate(alldata$alpha, alldata$Event)/alldata$Delta.time
#10. Subset data by side ####
pdr<-subset(alldata, Fish.Side=="RIGHT", row.names=NULL)
pdl<-subset(alldata, Fish.Side=="LEFT", row.names=NULL)
par(oldpar)
par(mfrow=c(3,1), mar=c(5,5,1,2))
plot(overallslope, type="l")
plot(alldata$x.filt, type="l")
plot(fs, type="l")
par(oldpar)
par(mfrow=c(1,1), mar=c(5,5,5,3))
plot(alldata$Time.diff, alldata$fish.temp, type="l", col="black", main="Behavioural thermoregulation",
xlab="Elapsed time (s)", ylab="Selected temperature (°C)")
points(pdr$Time.diff, pdr$fish.temp, col="red", cex=0.2)
points(pdl$Time.diff, pdl$fish.temp, col="blue", cex=0.2)
points(alldata$Time.diff[le.index], let, col="blue")
points(alldata$Time.diff[ue.index], uet, col="red")
par(mfrow=c(1,2), mar=c(5,5,5,5))
hist(let, main="Lower Escape Temperatures", col="blue")
hist(uet, main="Upper Escape Temperatures", col="red")
stat.desc(let)
## nbr.val nbr.null nbr.na min max
## 5.210e+02 0.000e+00 0.000e+00 9.400e+00 1.800e+01
## range sum median mean SE.mean
## 8.600e+00 6.916e+03 1.340e+01 1.327e+01 8.438e-02
## CI.mean.0.95 var std.dev coef.var
## 1.658e-01 3.709e+00 1.926e+00 1.451e-01
stat.desc(uet)
## nbr.val nbr.null nbr.na min max
## 5.210e+02 0.000e+00 0.000e+00 1.240e+01 2.110e+01
## range sum median mean SE.mean
## 8.700e+00 8.545e+03 1.650e+01 1.640e+01 8.837e-02
## CI.mean.0.95 var std.dev coef.var
## 1.736e-01 4.069e+00 2.017e+00 1.230e-01
mean(fish.temp)
## [1] 15.07
mean(let)
## [1] 13.27
mean(uet)
## [1] 16.4
par(mfcol=c(2,2), mar=c(5,5,2,2))
hist(alldata$velocity[le.index], col="blue")
hist(alldata$velocity[ue.index], col="red")
hist(alldata$acceleration[le.index], col="blue")
hist(alldata$acceleration[ue.index], col="red")
par(mfrow=c(1,1), mar=c(5,5,2,5))
smoothScatter(alldata$fish.temp, alldata$velocity, nbin=256,
colramp=heat.colors,
main="Swimming Velocity ~ Selected Temperature",
xlab="Selected Temperature",
"ylab"="Velocity",
postPlotHook=addgradient)
fit<-lm(velocity~fish.temp, data=alldata)
abline(coef=fit$coef, xpd=FALSE)
#11. Plot trajectory data based on side ####
par(mfrow=c(1,1), mar=c(5,5,2,2))
plot(pdr$Time.diff, pdr$velocity, type="p", col="red", main="Velocity by side",
xlab="Elapsed time (s)", ylab="Velocity (cm/s)", cex=0.2)
points(pdl$Time.diff, pdl$velocity, col="blue", cex=0.2)
plot(pdr$Time.diff, pdr$acceleration, type="p", col="red", main="Acceleration by side",
xlab="Elapsed time (s)", ylab="Acceleration (cm/s/s)", cex=0.2)
points(pdl$Time.diff, pdl$acceleration, col="blue", cex=0.2)
plot(pdr$Time.diff, pdr$alpha, type="p", col="red", main="Alpha by side",
xlab="Elapsed time (s)", ylab="Turning angle", cex=0.2)
points(pdl$Time.diff, pdl$alpha, col="blue", cex=0.2)
par(mfrow=c(1,1), xpd=TRUE)
polar.plot(pdr$velocity, pdr$theta, point.col="red", rp.type="s", start=0,
main="Movement vectors (° vs. cm/s)", cex=0.2, mar=c(4,4,4,4))
polar.plot(pdl$velocity, pdl$theta, point.col="blue", cex=0.2, rp.type="s", start=0,
add=TRUE)
legend(x="topright", c("Right side", "Left side"), col=c("red", "blue"), pch=c(1,1),
bty="n", inset=c(-0.2,0))
par(mfrow=c(1,1), xpd=TRUE)
polar.plot(pdr$acceleration[-1], pdr$alpha[-1], point.col="red", start=90, rp.type="s",
main="Turning vectors (° vs. cm/s/s)", mar=c(4,4,4,4), cex=0.2)
polar.plot(pdl$acceleration[-1], pdl$alpha[-1], point.col="blue", start=90, rp.type="s",
add=TRUE, cex=0.2)
legend(x="topright", c("Right side", "Left side"), col=c("red", "blue"), pch=c(1,1),
bty="n", inset=c(-0.2,0))
par(oldpar)
plot(cumsum(na.omit(alldata$alpha)), main="Cumulative Alpha")
plot(cumsum(na.omit(pdl$alpha)), main="Cumulative Alpha Left")
plot(cumsum(na.omit(pdr$alpha)), main="Cumulative Alpha Right")
plot(cumsum(na.omit(alldata$distance)), main="Cumulative Distance")
plot(cumsum(na.omit(pdl$distance)), main="Cumulative Distance Left")
plot(cumsum(na.omit(pdr$distance)), main="Cumulative Distance Right")
# if you take the cumulative sum of the instantaneous turning angles, you can derive an overall
# impression of whether one direction dominates more than the other.
# also, on each side of the shuttle box, the cumulative sum will indicate if turns occur more
# in one way or the other
library(ggplot2)
g<-ggplot(data=pdl, aes(x=1:nrow(pdl), y=cumalpha))+
geom_point(aes(col=alpha), alpha=1)+
scale_color_gradientn(colors=ironbowpal, limits=c(-50,50))
g
g<-ggplot(data=pdr, aes(x=1:nrow(pdr), y=cumalpha))+
geom_point(aes(col=alpha), alpha=1)+
scale_color_gradientn(colors=ironbowpal, limits=c(-50,50))
g
g<-ggplot(data=alldata, aes(x=x.filt, y=y.filt))+
geom_point(aes(col=alpha, size=abs(alpha)), alpha=1)+
scale_color_gradientn(colors=ironbowpal, limits=c(-50,50))
g
g<-ggplot(data=pdl, aes(x=x.filt, y=y.filt))+
geom_point(aes(col=alpha), size=1, alpha=0.5)+
scale_color_gradientn(colors=ironbowpal, limits=c(-20,20))
g
g<-ggplot(data=pdr, aes(x=x.filt, y=y.filt))+
geom_point(aes(col=alpha), size=1, alpha=0.5)+
scale_color_gradientn(colors=ironbowpal, limits=c(-20,20))
g
# Area Plot ####
# size of fish ~ Area
g<-ggplot(data=alldata, aes(x=x.filt, y=y.filt))+
ggtheme()+
geom_point(aes(col=Area, size=Area), alpha=0.2)+
scale_size_continuous(range=c(0.5,3))+
scale_color_gradientn(colors=ironbowpal)+
xlab("X (cm)")+ylab("Y (cm)")
g
#12. Summarise and Output #####
options(scipen=999)
colstats<-c("fish.temp", "Delta.time", "velocity", "theta", "acceleration", "alpha")
overall.sum<-stat.desc(alldata[colstats], basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
overall.sum
## fish.temp Delta.time velocity theta
## nbr.val 40969.00000 40969.0000000 40969.00000000000000 40969.0000
## nbr.null 0.00000 0.0000000 0.00000000000000 16.0000
## nbr.na 0.00000 0.0000000 0.00000000000000 0.0000
## min 9.40000 1.0000000 0.00000000005707 -180.0000
## max 21.10000 2.0000000 3.74689477794547 180.0000
## range 11.70000 1.0000000 3.74689477788840 360.0000
## sum 617610.00000 41616.0000000 25206.30479494316751 218597.2593
## median 15.10000 1.0000000 0.43742721507705 -2.0503
## mean 15.07506 1.0157924 0.61525311320616 5.3357
## SE.mean 0.01133 0.0006159 0.00305534675004 0.5242
## CI.mean.0.95 0.02221 0.0012073 0.00598854651678 1.0275
## var 5.25911 0.0155434 0.38245150482482 11258.5763
## std.dev 2.29328 0.1246732 0.61842663657448 106.1064
## coef.var 0.15212 0.1227349 1.00515807770850 19.8862
## acceleration alpha
## nbr.val 40968.00000000000000000 40968.0000
## nbr.null 0.00000000000000000 53.0000
## nbr.na 1.00000000000000000 1.0000
## min 0.00000000000005844 -180.0000
## max 1.62234725198670482 180.0000
## range 1.62234725198664642 360.0000
## sum 6790.79102449924630491 -36967.0621
## median 0.09643254418079467 -0.3733
## mean 0.16575842180480488 -0.9023
## SE.mean 0.00093462613372195 0.2148
## CI.mean.0.95 0.00183188768396192 0.4210
## var 0.03578661357096273 1889.8866
## std.dev 0.18917350123884352 43.4728
## coef.var 1.14126027009120512 -48.1779
right.sum<-stat.desc(pdr[colstats], basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
right.sum
## fish.temp Delta.time velocity theta
## nbr.val 21562.00000 21562.0000000 21562.00000000000000 21562.0000
## nbr.null 0.00000 0.0000000 0.00000000000000 16.0000
## nbr.na 0.00000 0.0000000 0.00000000000000 0.0000
## min 12.40000 1.0000000 0.00000000005707 -180.0000
## max 21.10000 2.0000000 3.70996714368215 180.0000
## range 8.70000 1.0000000 3.70996714362508 360.0000
## sum 351454.40000 21885.0000000 12453.46918036453280 37128.4729
## median 16.50000 1.0000000 0.42986567304395 -1.6639
## mean 16.29971 1.0149801 0.57756558669718 1.7219
## SE.mean 0.01385 0.0008273 0.00376830524324 0.7262
## CI.mean.0.95 0.02715 0.0016215 0.00738615719406 1.4233
## var 4.13574 0.0147563 0.30618308244764 11369.7492
## std.dev 2.03365 0.1214757 0.55333812668895 106.6290
## coef.var 0.12477 0.1196828 0.95805245228205 61.9238
## acceleration alpha
## nbr.val 21561.00000000000000000 21561.0000
## nbr.null 0.00000000000000000 53.0000
## nbr.na 1.00000000000000000 1.0000
## min 0.00000000000005844 -180.0000
## max 1.49025557583006996 180.0000
## range 1.49025557583001156 360.0000
## sum 3232.06811766507280481 3907.4290
## median 0.09026910555416823 0.0000
## mean 0.14990344221812871 0.1812
## SE.mean 0.00115862578906106 0.2879
## CI.mean.0.95 0.00227099231006296 0.5644
## var 0.02894378219702692 1787.4958
## std.dev 0.17012872243400562 42.2788
## coef.var 1.13492205326710605 233.2923
left.sum<-stat.desc(pdl[colstats], basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
left.sum
## fish.temp Delta.time velocity theta
## nbr.val 19407.00000 19407.0000000 19407.0000000 19407.0000
## nbr.null 0.00000 0.0000000 0.0000000 0.0000
## nbr.na 0.00000 0.0000000 0.0000000 0.0000
## min 9.40000 1.0000000 0.0002651 -180.0000
## max 18.10000 2.0000000 3.7468948 180.0000
## range 8.70000 1.0000000 3.7466296 360.0000
## sum 266155.60000 19731.0000000 12752.8356146 181468.7864
## median 13.70000 1.0000000 0.4465202 -2.5535
## mean 13.71441 1.0166950 0.6571256 9.3507
## SE.mean 0.01241 0.0009197 0.0048890 0.7564
## CI.mean.0.95 0.02433 0.0018028 0.0095829 1.4827
## var 2.98965 0.0164171 0.4638776 11105.0069
## std.dev 1.72906 0.1281293 0.6810856 105.3803
## coef.var 0.12608 0.1260254 1.0364619 11.2698
## acceleration alpha
## nbr.val 19407.000000000 19407.0000
## nbr.null 0.000000000 0.0000
## nbr.na 0.000000000 0.0000
## min 0.000007746 -180.0000
## max 1.622347252 180.0000
## range 1.622339506 360.0000
## sum 3558.722906834 -40874.4911
## median 0.106391432 -1.3378
## mean 0.183373160 -2.1062
## SE.mean 0.001485077 0.3211
## CI.mean.0.95 0.002910879 0.6294
## var 0.042801229 2000.9861
## std.dev 0.206884580 44.7324
## coef.var 1.128216258 -21.2387
alldata$TimeIndex<-floor(alldata$Time.diff/60)
tempsum<-ddply(alldata, .(TimeIndex), summarize, MeanTemp=mean(fish.temp, na.rm=TRUE))
velocitysum<-ddply(alldata, .(TimeIndex), summarize, MeanVelocity=mean(velocity, na.rm=TRUE))
thetasum<-ddply(alldata, .(TimeIndex), summarize, MeanTheta=mean(theta, na.rm=TRUE))
accelerationsum<-ddply(alldata, .(TimeIndex), summarize, MeanAcceleration=mean(acceleration, na.rm=TRUE))
alphasum<-ddply(alldata, .(TimeIndex), summarize, MeanAlpha=mean(alpha, na.rm=TRUE))
areasum<-ddply(alldata, .(TimeIndex), summarize, MeanArea=mean(Area, na.rm=TRUE))
# merge 3 data frames at once:
temp<-list(tempsum, velocitysum, thetasum, accelerationsum, alphasum, areasum)
fishtimesum<-Reduce(function(x, y) merge(x, y, all=TRUE, sort=F), temp)
# Finalise summary data
Expt.Duration<-NA
Time.low<-NA
Time.upp<-NA
Duration.low<-NA
Duration.upp<-NA
No.shuttles<-NA
LET.min<-NA
LET.max<-NA
LET.avg<-NA
LET.med<-NA
LET.sd<-NA
UET.min<-NA
UET.max<-NA
UET.avg<-NA
UET.med<-NA
UET.sd<-NA
Area.med<-NA
TFish.avg<-NA
Tfish.med<-NA
Tfish.sd<-NA
Vel.avg<-NA
Vel.med<-NA
Theta.avg<-NA
Theta.med <-NA
Accel.avg<-NA
Accel.med<-NA
Alpha.avg <-NA
Alpha.med<-NA
Vel.avg.f<-NA
Accel.avg.f<-NA
Theta.avg.f <-NA
Alpha.avg.f <-NA
CumulDist.Left <-NA
CumulDist.Right<-NA
Bias.index<-NA
Bias.strength<-NA
Expt.Duration<-max(alldata$Time.diff, na.rm=TRUE) # no. of seconds in the entire experiment
Time.low<-sum(pdl$Delta.time, na.rm=TRUE)
# no. of seconds spent in the left side
Time.upp<-sum(pdr$Delta.time, na.rm=TRUE)
# no. of seconds spent in the right side
Duration.low<-Time.low/nrow(leshuttles)
# average duration of time spent in the left side
Duration.upp<-Time.upp/nrow(ueshuttles)
# average duration of time spent in the right side
No.shuttles<-min(nrow(leshuttles),nrow(ueshuttles), na.rm=TRUE)
# number of shuttles overall. take the min value of lower and upper shuttles
LET.min<-min(leshuttles$LET, na.rm=TRUE)
# minumum lower escape temperature
LET.max<-max(leshuttles$LET, na.rm=TRUE)
# maximum lower escape temperature
LET.avg<-mean(leshuttles$LET, na.rm=TRUE)
# average lower escape temperature
LET.med<-median(leshuttles$LET, na.rm=TRUE)
# median lower escapte temperature (in case distribution is skewed)
LET.sd<-sd(leshuttles$LET, na.rm=TRUE)
UET.min<-min(ueshuttles$UET, na.rm=TRUE)
# minimum upper escape temperature
UET.max<-max(ueshuttles$UET, na.rm=TRUE)
# maximum upper escape temperature
UET.avg<-mean(ueshuttles$UET, na.rm=TRUE)
# average upper escape temperature
UET.med<-median(ueshuttles$UET, na.rm=TRUE)
# median upper escape temperature
UET.sd<-sd(ueshuttles$UET, na.rm=TRUE)
Area.med<-median(alldata$Area, na.rm=TRUE)
#### No filtering
TFish.avg<-mean(alldata$fish.temp, na.rm=TRUE)
# overall average temperature of the fish, based on total time spent at any given temperature
Tfish.med<-median(alldata$fish.temp, na.rm=TRUE)
# overall median temperature of the fish, based on total time spent at any given temperature
Tfish.sd<-sd(alldata$fish.temp, na.rm=TRUE)
Vel.avg<-mean(alldata$velocity, na.rm=TRUE)
# overall average velocity, includes all the 'stationary' periods
Vel.med<-median(alldata$velocity, na.rm=TRUE)
# overall median velocity, includes all the 'stationary' periods
Theta.avg<-mean(alldata$theta, na.rm=TRUE)
# overall average theta
Theta.med<-median(alldata$theta, na.rm=TRUE)
# overall median theta
Accel.avg<-mean(alldata$acceleration, na.rm=TRUE)
# overall average acceleration, includes all the 'stationary' periods
Accel.med<-median(alldata$acceleration, na.rm=TRUE)
# overall median acceleration, includes all the 'stationary' periods
Alpha.avg<-mean(alldata$alpha, na.rm=TRUE)
# overall average alpha, including the 'stationay' periods
Alpha.med<-median(alldata$alpha, na.rm=TRUE)
# overall median alpha, including the 'stationay' periods
CumulDist.Left<-sum(pdl$distance)
# cumulative distance moved when on the left side
CumulDist.Right<-sum(pdr$distance)
# cumulative distnace moved when on the right side
filtered.data<-subset(alldata, velocity>Vel.med & abs(alpha)>20)
# create new data according to filter of velocities>median and alphas>20 degrees
counterclockwise<-subset(filtered.data[colstats], alpha>20)
# subset only those alphas that are >20 degrees into counterclockwise dataframe
clockwise<-subset(filtered.data[colstats], alpha<(-20))
# subset only those alphas that are < -20 degrees into clockwise dataframe
n.counterclockwise<-nrow(counterclockwise)
n.clockwise<-nrow(clockwise)
Bias.index<-(n.clockwise-n.counterclockwise)/(n.clockwise+n.counterclockwise)
# bias.index is the ratio of the difference in turns to the total counted turns
# -1 = all left, counterclockwise turns
# +1 = all right, clockwise turns
# 0 = no difference in turns
Bias.strength<-abs(mean(clockwise$alpha, na.rm=T)+mean(counterclockwise$alpha, na.rm=T))
cat(paste("The bias index (-1 = counterclockwise, +1 = clockwise) is: \n",
round(Bias.index,4)), "\n")
## The bias index (-1 = counterclockwise, +1 = clockwise) is:
## -0.0486
#cat('\n')
cat(paste("The bias strength (|difference in mean ° counterclockwise vs. clockwise|) is: \n",
round(Bias.strength,4), "\n"))
## The bias strength (|difference in mean ° counterclockwise vs. clockwise|) is:
## 0.924
#cat('\n')
### With filtering
Vel.avg.f<-mean(filtered.data$velocity, na.rm=TRUE)
Accel.avg.f<-mean(filtered.data$acceleration, na.rm=TRUE)
Theta.avg.f<-mean(filtered.data$theta, na.rm=TRUE)
Alpha.avg.f<-mean(filtered.data$alpha, na.rm=TRUE)
Filename<-f
summaryoutput<-data.frame(Filename, Expt.Duration, Time.low, Time.upp, Duration.low,
Duration.upp, No.shuttles,
LET.min, LET.max, LET.avg, LET.med, LET.sd,
UET.min, UET.max, UET.avg, UET.med, UET.sd,
Area.med,
TFish.avg, Tfish.med, Tfish.sd, Vel.avg, Vel.med, Theta.avg, Theta.med,
Accel.avg, Accel.med, Alpha.avg, Alpha.med,
Vel.avg.f, Accel.avg.f, Theta.avg.f, Alpha.avg.f,
CumulDist.Left, CumulDist.Right,
Bias.index, Bias.strength)
t(summaryoutput)
## [,1]
## Filename "3118 - 15LSJ11.csv"
## Expt.Duration "41616"
## Time.low "19731"
## Time.upp "21885"
## Duration.low "37.87"
## Duration.upp "42.01"
## No.shuttles "521"
## LET.min "9.4"
## LET.max "18"
## LET.avg "13.27"
## LET.med "13.4"
## LET.sd "1.926"
## UET.min "12.4"
## UET.max "21.1"
## UET.avg "16.4"
## UET.med "16.5"
## UET.sd "2.017"
## Area.med "402"
## TFish.avg "15.08"
## Tfish.med "15.1"
## Tfish.sd "2.293"
## Vel.avg "0.6153"
## Vel.med "0.4374"
## Theta.avg "5.336"
## Theta.med "-2.05"
## Accel.avg "0.1658"
## Accel.med "0.09643"
## Alpha.avg "-0.9023"
## Alpha.med "-0.3733"
## Vel.avg.f "0.7978"
## Accel.avg.f "0.2592"
## Theta.avg.f "1.635"
## Alpha.avg.f "1.529"
## CumulDist.Left "12857"
## CumulDist.Right "12538"
## Bias.index "-0.04862"
## Bias.strength "0.924"
setwd(outputDir)
if(file.exists(fileout))
{
write.table(summaryoutput, fileout, append=TRUE, sep=",", col.names=FALSE, row.names=FALSE)
}
if(!file.exists(fileout))
{
write.table(summaryoutput, fileout, append=TRUE, sep=",", col.names=TRUE, row.names=FALSE)
}
# Final File prep with processed data
processed<-list(overall.sum, summaryoutput, pdl, pdr, filtered.data, clockwise,
counterclockwise, leshuttles, ueshuttles, allshuttles, fishtimesum)
names(processed)<-c("overall.sum", "summaryoutput", "pdl", "pdr", "filtered.data", "clockwise",
"counterclockwise", "leshuttles", "ueshuttles", "allshuttles", "fishtimesum")
saveRDS(processed, file=Rdatafileout)
p<-readRDS(Rdatafileout)
write.csv(fishtimesum, file=timesummaryout, row.names=F )
write.csv(allshuttles, file=shuttlesummaryout, row.names=F)
#} # this closes the m loop used for randomisation procedure
#} # this closes the f in l.files[] loop